{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "# 寻找线性关系\n",
    "def Search_linear_relations(y, r, t):\n",
    "    time_start = time.clock()\n",
    "    L = Matrix(ZZ, r, r + t)\n",
    "    C = Matrix(ZZ, 3, r)\n",
    "    for i in range(0, r):\n",
    "        for j in range(0, r + t):\n",
    "            L[i, j] = 0\n",
    "        L[i, i + t] = 1\n",
    "        for j in range(0, t):\n",
    "            L[i, j] = y[i + j]       \n",
    "    reduced_L = L.BKZ(algorithm = 'NTL', fp = 'fp', block_size = 20)\n",
    "    \n",
    "    for i in range(0, 3):\n",
    "        for j in range(0, r):\n",
    "            C[i, j] = reduced_L[i, j + t]\n",
    "    \n",
    "    time_end = time.clock()\n",
    "    print('Step 1: time of searching polynomials that annihilate the sequence a: ', time_end - time_start, 's\\n\\n')\n",
    "    return C\n",
    "\n",
    "# 结式法还原模数\n",
    "def Recover_the_modulus(C, r):\n",
    "    time_start = time.clock()\n",
    "    R.<x> = ZZ[]\n",
    "    f1 = C[0, 0]\n",
    "    f2 = C[1, 0]\n",
    "    f3 = C[2, 0]\n",
    "    for i in range(1, r):\n",
    "        f1 = f1 + x^i * C[0, i]\n",
    "        f2 = f2 + x^i * C[1, i]\n",
    "        f3 = f3 + x^i * C[2, i]\n",
    "    \n",
    "    r1 = f1.resultant(f2)\n",
    "    r2 = f2.resultant(f3)\n",
    "    r3 = f3.resultant(f1)\n",
    "    \n",
    "    m_power = gcd(gcd(r1, r2), r3)\n",
    "    m = m_power^(1/16)\n",
    "    \n",
    "    time_end = time.clock()\n",
    "    print('Step 2: time of recovering the modulus: ', time_end - time_start, 's\\n\\n')\n",
    "    return m\n",
    "\n",
    "# GCD法还原系数\n",
    "def Recover_the_coefficients(C, m, r):\n",
    "    time_start = time.clock()\n",
    "    R.<x> = PolynomialRing(Integers(m), implementation='NTL')\n",
    "    \n",
    "    for i in range(0, 2):\n",
    "        mark = r - 1\n",
    "        while C[i][mark] == 0:\n",
    "            mark = mark - 1\n",
    "        for j in range(0, r):\n",
    "            C[i, j] = C[i, j] * ((C[i, mark]^(-1)) % m)\n",
    "    f = C[0, 0]\n",
    "    g = C[1, 0]\n",
    "    for i in range(1, r):\n",
    "        f = f + x^i * C[0, i]\n",
    "        g = g + x^i * C[1, i]\n",
    "        \n",
    "    time_end = time.clock()\n",
    "    print('Step 3: time of recovering the coefficients: ', time_end - time_start, 's\\n\\n')\n",
    "    return f.gcd(g)\n",
    "\n",
    "# Kannan嵌入法还原初态\n",
    "def Recover_the_differences_of_initial_state(y, f, m, n, beta):\n",
    "    time_start = time.clock()\n",
    "    d = 30\n",
    "    z = vector(ZZ, n)\n",
    "    a = vector(ZZ, n)\n",
    "    Q = Matrix(ZZ, n, n)\n",
    "    Q_power = identity_matrix(ZZ, n)\n",
    "    L = Matrix(ZZ, d + 1, d + 1)\n",
    "    \n",
    "    Q[0, n - 1] = (-f.list()[0]) % m\n",
    "    for i in range(1, n):\n",
    "        Q[i, i - 1] = 1\n",
    "        Q[i, n - 1] = (-f.list()[i]) % m\n",
    "    \n",
    "    for i in range(1, n):\n",
    "        Q_power = (Q_power * Q) % m\n",
    "    for i in range(n, d):\n",
    "        Q_power = (Q_power * Q) % m\n",
    "        b = 0\n",
    "        for j in range(0, n):\n",
    "            L[j + 1, i + 1] = Q_power[j, 0]\n",
    "            b = b + Q_power[j, 0] * y[j]\n",
    "        L[0, i + 1] = (2^beta * (y[i] - b)) % m + 2^(beta - 1)\n",
    "    \n",
    "    L[0, 0] = 2^(beta - 1)\n",
    "    for i in range(1, n + 1):\n",
    "        L[0, i] = 2^(beta - 1)\n",
    "        L[i, i] = 1\n",
    "    for i in range(n + 1, d + 1):\n",
    "        L[i, i] = m\n",
    "        \n",
    "    reduced_L = L.BKZ(algorithm = 'NTL', fp = 'fp', block_size = 20)\n",
    "    \n",
    "    if reduced_L[0, 0] == -2^(beta - 1):\n",
    "        for i in range(0, n):\n",
    "            z[i] = reduced_L[0, i + 1] + 2^(beta - 1)\n",
    "            a[i] = y[i] * 2^(beta) + z[i]\n",
    "    if reduced_L[0, 0] == 2^(beta - 1):\n",
    "        for i in range(0, n):\n",
    "            z[i] = -reduced_L[0, i + 1] + 2^(beta - 1)\n",
    "            a[i] = y[i] * 2^(beta) + z[i]\n",
    "            \n",
    "    time_end = time.clock()\n",
    "    print('Step 4: time of recovering the initial state: ', time_end - time_start, 's\\n\\n')\n",
    "    return a\n",
    "    \n",
    "    \n",
    "\n",
    "n = 16\n",
    "r = 175\n",
    "t = 65\n",
    "N = r + t - 1\n",
    "alpha = 17\n",
    "beta = 14\n",
    "\n",
    "y = vector([312,6195,20905,36760,33275,128773,5120,82629,3434,13450,115969,78595,46080,21032,126472,116368,39147,118169,53608,16848,30654,130814,92427,66994,64940,95674,83420,12319,51686,121628,103101,4578,55561,114165,99105,66029,11183,33637,126822,45172,71211,75549,31249,100570,57872,40625,53029,98594,35382,114241,86331,100286,74451,71452,113487,22824,119007,42576,9420,18312,62834,3099,103626,124788,28318,119993,120470,1417,66987,121877,100551,36538,107257,95791,90747,11856,78683,27967,33439,112227,42492,76803,73986,127640,127988,113849,62073,121003,67257,62990,28618,14066,37196,88819,91188,119086,47509,14040,46365,117363,96302,61501,69537,19261,5302,29632,112053,60563,73212,123325,7485,52451,92665,108139,31895,6630,28136,97771,1874,10433,123318,98757,45897,34866,118234,75543,95739,41629,82990,25310,28330,21183,102442,39164,13744,125171,78275,42927,27280,113230,115192,108433,91178,101247,100622,63782,105110,40546,20787,97620,33627,97857,54439,74890,65682,117061,34213,28237,72077,40838,44105,41674,48232,10965,79577,106188,27672,79314,114402,64028,44159,130067,47639,20721,16188,1301,53216,29228,52467,45721,17559,28506,122718,84634,109576,34973,129964,10245,48392,100922,5358,128436,37415,74936,89446,101529,75414,39696,7480,20690,113561,65658,28894,48013,100734,80683,119101,50176,38679,50948,30958,67432,120501,62535,103920,69944,73629,7114,19565,10250,685,37711,17024,38274,53345,89835,93532,61993,35816,55362,11054,112389,34513,44560,123360,26732,69001,67919,67745])\n",
    "C = Search_linear_relations(y, r, t)\n",
    "m = Recover_the_modulus(C, r)\n",
    "print('m = ',m)\n",
    "f = Recover_the_coefficients(C, m, r)\n",
    "print('f = ',f)\n",
    "a = Recover_the_differences_of_initial_state(y, f, m, n, beta)\n",
    "print('a = ',a)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.1",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
